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ABSTRACT 

We present the results of global cylindrical disc simulations and local shearing box 
simulations of protoplanets interacting with a disc undergoing MHD turbulence. The 
specific emphasis of this paper is to examine and quantify the magnitude of the torque 
exerted by the disc on the embedded protoplanets as a function of the protoplanet 
mass, and thus to make a first study of the induced orbital migration of protoplanets 
resulting from their interaction with magnetic, turbulent discs. This issue is of crucial 
importance in understanding the formation of gas giant planets through the so-called 
core instability model, and the subsequent orbital evolution post formation prior to 
the dispersal of the protostellar disc. Current estimates of the migration time of proto- 
planetary cores in the 3-30 Earth mass range in standard disc models are Tmig — 10* 
- 10^ yr, which is much shorter than the estimated gas accretion time scale of Jupiter 
type planets. 

The global simulations were carried out for a disc with aspect ratio H/R ~ 0.07 
and protoplanet masses of Mp = 3, 10, 30 Earth masses, and 3 Jupiter masses. The 
local shearing box simulations were carried out for values of the dimensionless param- 
eter {Mp/M^)/{H/Rf = 0.1, 0.3, 1.0, and 2.0, with being the central mass. These 
allow both embedded and gap forming protoplanets for which the disc response is non 
linear to be investigated. 

In all cases the instantaneous net torque experienced by a protoplanet showed 
strong fluctuations on an orbital time scale, and in the low mass embedded cases 
oscillated between negative and positive values. Consequently, in contrast to the lam- 
inar disc type I migration scenario, orbital migration would occur as a random walk. 
Running time averages for embedded protoplanets over typical run times of 20 - 25 
orbital periods, indicated that the averaged torques from the inner and outer disc took 
on values characteristic of type I migration. However, large fluctuations occurring on 
longer than orbital time scales remained, preventing convergence of the average torque 
to well defined values or even to a well defined sign for these lower mass cases. 

Fluctuations became relatively smaller for larger masses indicating better conver- 
gence properties, to the extent that in the 3OM0 simulation consistently inward, albeit 
noisy, migration was indicated. 

Both the global and local simulations showed this trend with increasing proto- 
planet mass which is due to its perturbation on the disc increasing to become compa- 
rable to and then dominate the turbulence in its neighbourhood. This then becomes 
unable to produce very large long term fluctuations in the torques acting on the pro- 
toplanet. Eventually gap formation occurs and there is a transition to the usual type 
H migration at a rate determined by the angular momentum transport in the distant 
parts of the disc. 

The existence of significant fiuctuations occurring in the turbulent discs on long 
time scales is an important unexplored issue for the lower mass embedded protoplanets, 
that are unable to modify the turbulence in their neighbourhood, and which have been 
studied here. If significant fiuctuations occur on the longest disc evolutionary time 
scales, convergence of torque running averages for practical purposes will not occur 
and the migrational behaviour of low mass protoplanets considered as an ensemble 
would be very different from predictions of type I migration theory for laminar discs. 
The fact that noise levels were relatively smaller in the local simulations may indicate 
the presence of long term global fluctuations, but the issue remains an important one 
for future investigation. 
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1 INTRODUCTION 

The ongoing discovery of extrasolar giant planets has 

stimulated renewed interest in the theory of planet forma- 
tion (e.g. Mayor & Queloz 1995; Marcy, Cochran, & Mayor 
1999; Vogt ct al. 2002; Santos et al. 2003). In the most com- 
monly accepted theory of how planets form, the so-called 
core instability model, gas giant planets form through the 
build-up of a rocky and icy core of ^ 15 Earth masses, which 
then undergoes gas accretion resulting in a gas-giant planet 
(e.g. Bodenheimer & Pollack 1986; Pollack et al. 1996). An 
alternative model involves the formation of giant planets 
through the gravitational fragmentation of the protostellar 
disc during the earlier phases of its evolution (Boss 2001). 
In either case disc-planet interaction will play an important 
role in the subsequent evolution. 

The gravitational interaction between protostellar discs 
and embedded protoplaricts has been the subject of a large 
number of studies over the last couple of decades. In the 
standard picture, a protoplanet exerts torques on a proto- 
stellar disc through the excitation of spiral density waves 
at Lindblad resonances, and possibly through interaction at 
corotation resonance (e.g. Goldroich & Tremaine 1979; Lin & 
Papaloizou 1979; Papaloizou & Lin 1984; Ward 1986, 1997; 
Tanaka, Tacheuchi & Ward 2002). The spiral waves carry 
with them an associated angular momentum flux. This an- 
gular momentum is deposited in the disc material where the 
waves are damped, leading to an exchange of angular mo- 
mentum between protoplanet and disc. The disc that lies 
exterior to the protoplanet orbit exerts a negative torque 
on the planet, and the interior disc exerts a positive torque. 
For most disc models the negative torque dominates and 
the protoplanet migrates inwards. For protoplanets of ~ 15 
Earth masses, the migration time is estimated to be between 
10" and 10^ yr (Tanaka, Tacheuchi, & Ward 2002), which is 
very much shorter than the estimated gas accretion phase of 
giant planet formation ~ 7 Myr (Pollack et al. 1996). Taken 
at face value, this presents a serious problem for the core- 
instability model of gas giant planet formation. We note, 
however, this analysis pertains only to smooth, laminar disc 
models. 

For protoplanets in the .Jovian mass range, the interac- 
tion is non linear and gap formation occurs (Papaloizou & 
Lin 1984; Bryden et al. 1999; Kley 1999). In this case the or- 
bital migration of the planet becomes locked to the viscous 
evolution of the disc, and migration is expected to occur on 
a time scale of 10^ yr (Lin & Papaloizou 1986; Nelson et al. 
2000; D'Angelo, Kley & Henning 2002). 

Until quite recently most models of viscous accretion 
discs used the Shakura & Sunyaev (1973) a model for the 
anomalous disc viscosity. This assumes that the viscous 
stress is proportional to the thermal pressure in the disc, 
without specifying the origin of the viscous stress (but as- 
sumed to arise from some form of turbulence). Work by 
Balbus & Hawley (1991) indicated that significant angu- 
lar momentum transport in weaJily magnetised discs could 
arise from the magnetorotational instability (MRI - or the 
Balbus-Hawlcy instability). Subsequent non linear numeri- 
cal simulations performed using a local shearing box formal- 
ism (e.g. Hawley & Balbus 1991; Hawley, Gammie, & Balbus 
1996; Brandenburg et al. 1996) confirmed this and showed 
that the saturated non linear outcome of the MRI is MHD 



turbulence with an associated viscous stress parameter a of 
between ~ 5 x 10~^ and ~ 0.1, depending on the initial 
magnetic field configuration. More recent global simulations 
of MHD turbulent discs [e.g. Armitage 1998; Hawley 2000; 
Hawley 2001; Steinacker & Papaloizou 2002; Papaloizou & 
Nelson 2003] confirm the picture provided by the local shear- 
ing box simulations. 

This is the fourth in a series of papers that examine the 
interaction between disc models undergoing MHD turbu- 
lence with zero net fiux magnetic fields and embedded pro- 
toplanets. In Papaloizou & Nelson (2003 - hereafter paper I) 
wo examined and characterised the turbulence obtained in a 
variety of MHD disc models. In Nelson & Papaloizou (2003 
- hereafter paper II) we examined the interaction between a 
global cylindrical disc model and a massive (5 Jupiter mass) 
protoplanet. A similar study was undertaken by Winters, 
Balbus, & Hawley (2003b). In a companion paper to this 
one (Papaloizou, Nelson, & Snellgrove 2003 - hereafter pa- 
per HI) we presented the results of global cylindrical disc 
simulations and local shearing box simulations of turbulent 
discs interacting with protoplanets of different mass. The 
main focus of paper III was to characterise the changes in 
fiow morphology and disc structure as a function of planet 
mass, and to examine the transition from linear to non lin- 
ear interaction leading to gap formation. In this paper we 
continue to examine these simulations, but now focus on the 
gravitational torques exerted on the protoplanet by the disc 
and the associated migration rate of the protoplanet. 

We find that in all simulations performed, the torque 
experienced by the protoplanet is a highly variable quantity 
on account of the protoplanet interacting with the turbulent 
density wakes that shear past it. For low mass protoplanets, 
the torque is dominated by these fluctuations, such that the 
usual distinction between inner (positive) and outer (nega- 
tive) disc torques is blurred. The net torcjuc experienced by 
embedded protoplanets oscillates between negative and pos- 
itive values, such that the protoplanet migration is likely to 
occur as a random walk. This is in contrast to the monotonic 
inward drift normally associated with type I migration. A 
running time average of the torques fails to converge for the 
embedded protoplanet runs, at least for the run times that 
are currently feasible, so that definitive statements about 
the direction and rate of migration of low mass planets in 
turbulent discs cannot yet be made. 

In a manner that is consistent with the results of paper 
HI, we find that the results show a definite trend as a func- 
tion of planet mass. For very low mass planets the turbulent 
density wakes are of much higher amplitude than the spiral 
wakes generated by the planet. This is reflected in the torque 
experience by the planet, for which the turbulent fluctua- 
tions are very much larger than the running mean torque. 
As the planet mass increases, the spiral wakes generated by 
the planet become apparent in the flow. This is accompanied 
by the torques displaying the expected separation between 
inner and outer disc torques, and the fluctuations becom- 
ing smaller relative to the running mean torque. This trend 
with increasing mass continues until gap formation occurs. 
At this point wo find a weakening of the torques exerted by 
the disc on the protoplanet due to the evacuation of gap 
material. There is then a transition to type II migration 
and, as exemplified by the simulation of a 5 Jupiter mass 
protoplanet interacting with a turbulent disc presented in 
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paper II, inward migration at a rate determined by the an- 
gular momentum transport in the distant parts of the disc 
unaffected by the protoplanet. 

The plan of the paper is as follows. In section 2 wc 
described our initial model set up and numerical procedure. 
In section 3 we present the results of the simulations. Finally 
in section 4 we discuss our results. 



2 INITIAL MODEL SETUP 

The initial set up and boundary conditions used in the mod- 
els presented in this paper are described in detail in paper 
III, and wc do not repeat that discussion hero. However, we 
present some details of the runs in tables 1 and 2 for con- 
venience. Table 1 contains details of the global runs, and 
table 2 contains details of the local shearing box runs. 

2.1 Global Runs 

The initial conditions for the global runs consisted of a tur- 
bulent accretion disc model with H/r — 0.07, H being the 
putative disc semi thickness, and a volume averaged char- 
acteristic Shakura-Sunyaev a value of a ~ 7 x 10~^, as 
described in paper III. At time t = the protoplanet was 
inserted into the disc and subsequently held fixed at loca- 
tion {rp,4>p) = (3, tt) in cylindrical coordinates (r,(f)) cen- 
tred on the central mass M, . The gravitational potential 
of the protoplanet was softened using a softening parameter 
b = 0.3Hp, where Hp is the disc scale height at the planet lo- 
cation. The simulation was performed in a rotating reference 
frame whose angular velocity equalled the Koplorian angular 
velocity of the protoplanet. The torque on the planet as a 
function of time was calculated in the manner described in 
section 2.4. The unit of time used when discussing the global 
runs is the orbital period of the protoplanet Pp = 2ir/flp. 



code. The code has been developed from a version of NIR- 
VANA originally written by U. Ziegler (Ziegler & Rudiger 
2000). 

2.4 Torque Calculation 

The torque oxporicncod by the protoplanet in the global disc 
runs was calculated by summing the gravitational force due 
to the mass in each grid cell in the active domain of the 
disc model. Those cells that lay within the boundary layers 
located near the radial boundaries of the computational do- 
main (described in paper III) did not contribute to the gravi- 
tational force. The gravitational force due to the protoplanet 
acting on the disc was softened using b = 0.3Hp where Hp is 
the putative scale height at the position of the protoplanet. 
An identical softening was used when calculating the force of 
the disc on the protoplanet. Material that lay inside the Hill 
sphere of the protoplanet given by Rh = rp{Mp/iMtY^^ 
was excluded from the torque calculation. The acceleration 
experienced by the protoplanet can be written as 

J _ G nii (vp - Vj) 

where the sum is over all possible values of the subscript 
i which is used to label quantities evaluated at the centre 
of a specific grid cell, thus ensuring that the whole of the 
computational domain is covered. The mass in each grid 
cell rrii = pi Ar n A(f> Az, where Ar, A(f>, and Az are the 
grid spacings in the radial, azimuthal, and vertical directions 
respectively. The resulting torque per unit mass is then T = 
Vp X f . This quantity was calculated and stored every 10 time 
steps. 

A similar procedure was used to calculate the force act- 
ing on the planet in the shearing box simulations. 



2.2 Shearing Box Runs 

The morphology of the shearing box runs has already been 
described in detail in paper HI. These were generated from 
a simulation BaO described there which had fully developed 
turbulence, but no protoplanet, that had been run for 354 
time units. The time unit was taken to be the inverse angular 
velocity ,Qp, at the centre of the box located a putative dis- 
tance R from the central object. Simulations Bal-Ba4 were 
then carried out after inserting protoplanets with a range of 
values for GMp/{H^Ql) = MpR^/{M,H^), with Mp being 
the protoplanet mass. These and other relevant parameters 
are given in table 2. As for the global simulations, the pro- 
toplanet potential was softened, by use of a softening pa- 
rameter b = 0.3H. For additional details see paper HI. The 
results of these simulations are presented in section 3.6. 

2.3 Numerical procedure 

The numerical scheme that we employ is based on a spa- 
tially second-order accurate method that computes the ad- 
vection using the monotonic transport algorithm (Van Leer 
1977). The MHD section of the code uses the method of 
characteristics constrained transport (MOCCT) as outlined 
in Hawley & Stone (1995) and implemented in the ZEUS 



3 NUMERICAL RESULTS 

Before discussing the results of the global turbulent disc sim- 
ulations, we present the results of some simulations designed 
to calibrate the code. In particular, wc will examine the role 
that gravitational softening of the protoplanet gravitational 
potential has on the migration time in laminar disc simular- 
tions. 

3.1 Type I migration in laminar discs 

The migration time of a non gap forming protoplanet em- 
bedded in a three dimensional disc has been calculated using 
a linear analysis by Tanaka, Takeuchi, & Ward (2002) . They 
derive the following expression for the migration time of a 
protoplanet embedded in a locally isothermal disc with a 
power law surface density profile: 

Here Ep is the surface density at the position of the pro- 
toplanet, fip is the Keplerian angular velocity of the proto- 
planet, c is the sound speed in the disc at the protoplanet 
position, and —7 is the power law index for the radial surface 
density distribution. These authors comment that in general 
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Model 


cj) domain 


H/r 




(MpR^/{M,.H^) 
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n,f, 




Gl 
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1 X 10""' 


0.03 
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3 X lO-"' 


0.09 
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1092 


40 


G3 
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0.07 


1 X 10""' 


0.30 
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1092 


40 
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0.07 


3 X 10-3 


8.75 


450 


1092 


40 


G5 


7r/2 


0.07 


3 X 10-3 


8.75 
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276 


40 



Table 1. Parameters of the global simulations: The first column gives the simulation label, the second gives the extent of the aaimuthal 
domain, the third gives the H/r value of the disc, and the fourth gives the protoplanet-star mas ratio. The fifth gives the ratio of Mp/M« 
to {H/rY . The sixth, seventh, and eighth columns describe the number of grid cells used in each coordinate direction. 
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Table 2. Parameters of the shearing box simulations: The first column gives the simulation label, the second, third and fourth give the 
extent of the coordinate domains considered. The x, y, and z domains refer to the Cartesian domains labeled as radial, azimuthal and 
vertical respectively. The fifth and sixth column give the start and end times if the simulation measured in dimensionless units described 
in section 2.2. These were all started by inserting the protoplanet into the turbulent model generated from the simulation BaO described 
in paper III at time indicated. The seventh column gives the value of {GMp/(illH^) = {MpR^ /(M^H^). The eighth, ninth and tenth 
columns give the number of computational grid points in the x, y, and z coordinates respectively. These numbers include any ghost zones 
used to handle boundary conditions. 



this expression gives a migration time that is a factor of be- 
tween 2-3 times slower than similar expressions derived for 
flat, two-dimensional discs (e.g. Ward 1997). The finite disc 
thickness thus acts to soften the disc-planet gravitational 
interaction. 

The global turbulent disc simulations presented here 
use cylindrical disc models in which the vertical component 
of gravity is neglected, and thus the full three dimensional 
structure of the disc is not modeled. In the limit of a laminar 
disc, the cylindrical disc can be viewed as a series of two di- 
mensional discs stacked on top of one another, providing the 
protoplanet gravitational potential is also cylindrical (i.e. 
no vertical component). Given that gravitational softening 
of the protoplanet potential is employed, we are interested 
in (i) How large a gravitational softening is required in a 
two dimensional disc simulation to obtain agreement with 
equation 1; (ii) What is the deviation from equation 1 ob- 
tained for differing softening parameters. In order to address 
these points a number of two dimensional simulations using 
laminar, viscous a-discs have been performed with varying 
physical and softening parameters. 

Figure 1 shows the migration time for a series of models 
in which protoplanots of differing masses were placed in a 
protostellar disc model and the torque of the disc acting 
on the protoplanet was calculated. These torques were then 
used to calculate a migration time using the expression 



1 jp 
2T„ 



(2) 
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Figure 1. This figure shows the migration time for protoplan- 
ets of different mass using a gravitational softening parameter 
b = 0.7Hp. The dots show the results of the simulations, and the 
dashed line shows the migration time computed from equation 1. 
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Figure 2. This figure shows the migration time for protoplan- 
ets of different mass using a gravitational softening parameter 
b = 0.3Hp. The dots show the results of the simulations, and the 
dashed line shows the migration time computed from equation 1. 



where jp is the specific angular momentum of the proto- 
planet and Tp is the torque per unit mass due to the disc. 
The disc model used for the particular set of simulations 
presented in figure 1 was such that the disc surface den- 
sity E(r) — T,or~^^^, H/r = 0.05, and the dimensionless 
viscosity coefficient a = 4 x 10""^. In code units the inner 
boundary was located at Rin = 1, the outer boundary at 
Rout ~ 8, and the protoplanet was located at = 3. If 
we use the convention that Vp is equivalent to 5.2 AU, and 
the central stellar mass is a solar mass, then the disc mass 
was normalised so that it contained the equivalent of 7.5 
Jupiter masses between 1.56 and 20.8 AU. This disc model 
is equivalent to that described in Bate et al. (2003), and fig- 
ure 1 is directly comparable with their figure 10. The gravi- 
tational softening parameter was set to be 6 = Q.lHp where 
Hp is the disc thickness at the position of the protoplanet. 
We use the convention that 1 Earth mass corresponds to 
Mp/M. = 3 X 10"®. 

The black dots in figure 1 show the migration time 
scales obtained in the simulation. The dashed line is ob- 
tained from equation 1 assuming an identical disc model. It 
is clear that a gravitational softening parameter of 6 = 0.7-ffp 
gives good agreement with the migration time appropriate 
to a fully three dimensional locally isothermal disc, and that 
two dimensional simulations can be used to study the mi- 
gration of low mass protoplanets provided an appropriate 
gravitational softening parameter is used. 

Figure 2 shows the migration time calculated for simu- 
lations in which E(r) oc , H/r — 0.07, and the softening 
parameter was b — 0.3Hp. These are the values adopted for 
the global turbulent disc runs presented in sections 3.2, 3.3, 
and 3.4, because the potential is better represented close to 
the planet. The disc model here was normalised so that it 
contained the equivalent of 0.02 M0 between and 40 AU, 



Figure 3. This figure shows midplane density contours for the 
run Gl. Note that the presence of the protoplanet is undetectable 
due to the higher amplitude perturbations generated by the tur- 
bulence. The protoplanet is located at {xp,yp)= (-3,0). 

and ~ 2.5 Jupiter masses interior to the protoplanet ra- 
dius Vp. Figure 2 shows that this results in migration times 
that are a factor of ~ 2/3 shorter than those predicted by 
Tanaka, Tacheuchi, & Ward (2002) for the smaller masses 
but are in slightly closer agreement in the higher ~ 3OM0 
range. These shorter migration times arise because of the 
smaller softening parameter employed. 



3.2 Global Model Gl 

The global model Gl described in table 1 contained a proto- 
planet whose mass is equivalent to ~ 3 Earth masses. Such 
a protoplanet is expected to provide a small perturbation 
to the disc. A snapshot of the midplane density is shown 
in figure 3, and shows that the presence of the protoplanet 
is indiscernible due to the higher amplitude density fluctua- 
tions generated by the turbulence. Figure 4 shows a series of 
close-up images of the midplane density in the near vicinity 
of the protoplanet, which again show that the protoplanet 
cannot be detected. A similar plot for an equivalent laminar 
disc run is presented in figure 5 for purposes of comparison. 
It is worth noting that once the density structure shown in 
figure 5 is established in the laminar disc runs, it remains es- 
sentially time independent. The changing density structure 
observed in figure 4, combined with the larger amplitude of 
the turbulent wakes when compared with those generated 
by the protoplanet, suggest that the gravitational interac- 
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Mp=3 Earth Masses Mp=3 Earth Masses 




Figure 4. This figure shows midplane density contours for the run Gl in the region close to the protoplanet. The protoplanet is located 
at {rp,flip)= (3,7r)- The panels correspond to times 452.4903, 497.7334, 542.9972, and 672.4058, respectively. Note that the presence of 
the protoplanet is undetectable due to the higher amplitude perturbations generated by the turbulence. 
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Laminar Disc - Mp=3 Eorth Masses 
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Figure 5. This figure shows midplane density contours for a lami- 
nar disc run with planet parameters identical to run Gl. The pres- 
ence of the planet due to the excitation of spiral waves is clearly 
detectable, in contrast to the picture presented in figure 4. 



tion between a low mass protoplanet and turbulent disc may 
differ substantially when compared with a laminar disc. 

Figure 6 shows the time evolution of the torque per unit 
mass exerted on a protoplanet of equivalent parameters to 
that in run Gl but embedded in a laminar disc. The upper 
line shows the (positive) torque due to the disc lying interior 
to the protoplanet orbital radius, the lowest line shows the 
torque due to the outer disc, and the middle line shows the 
total (negative) torque. It is clear that a well defined net 
torque acting on the protoplanet may be defined, along with 
a corresponding inward migration rate. 

Figure 7 shows a similar plot but for the turbulent disc 
model Gl. It is clear that the forcing experienced by the pro- 
toplanet in this case differs dramatically from that obtained 
in the laminar disc. The torques from the inner and outer 
disc suffer very large fluctuations as a result of the proto- 
planet interacting gravitationally with the turbulent wakes 
apparent in figures 3 and 4 as they shear past the proto- 
planet, and this causes the net torque experienced by the 
protoplanet to oscillate between negative and positive val- 
ues. The orbital migration of the protoplanet is thus likely 
to occur as a random walk rather than as a monotonic in- 
ward drift normally associated with type I migration. In 
this run, for which the spiral wakes generated by the planet 
are completely dominated by the turbulent wakes, the usual 
separation between inner and outer disc torques is barely 
discernible due to the high amplitude turbulent fiuctuations. 



5 10 15 20 

Time (Orbits) 

Figure 6. This figure shows the torque per unit mass exerted by 
the disc on the protoplanet for a laminar disc and planet param- 
eters identical to run Gl. The upper line shows the torque due to 
the inner disc, the lower line shows the torque due to the outer 
disc, and the middle line shows the total torque. It is clear that a 
well defined net torque is produced, with an associated migration 
time scale. 



Turbulent Torques Per Unit Moss Vs Time (Mp = 3) 




Time (Orbits) 

Figure 7. This figure shows the torque per unit mass exerted by 
the disc on the protoplanet for the run Gl. The turbulence in this 
case generates strong fluctuations in the torque from each side of 
the disc, such that it becomes difficult to distinguish the torques 
arising from each side. The total torque fluctuates between pos- 
itive and negative values such that the associated migration will 
undergo a 'random walk'. 

The time evolution of the running time average of the 
torque per unit mass in model Gl is shown in figure 8. The 
upper line shows the running time average of the torque due 
to the inner disc, the lowest line shows the time averaged 
torque due to the outer disc, and the middle line shows the 
running time average of the total torque. The straight line 
shows the time average of the total torque per unit mass 
experienced in an equivalent laminar disc model. The first 
thing to note from figure 8 is that the running time averaged 
torque in run Gl does not converge to a well defined value 
for the time over which the simulation was run (just over 
20 planet orbits). Near the beginning of the simulation the 
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Time Averaged Torques Per Unit Mass (Mp = 3) 
-6 [ 




Figure 8. This figure shows the running time average of the 
torque per unit mass exerted by the disc for run Gl. The upper 
line is the running time average of the torque acting on the planet 
due to the inner disc. The lower line is that due to the outer disc. 
The middle (not straight) line is the running time average of the 
total torque. The straight line is the total torque exerted on the 
protoplanet in a laminar disc run. We note that the total time 
averaged torque does not converge to a well defined value. 



time average is close to that obtained from a similar laminar 
disc run, but as the calculation evolves the running time 
average changes continuously. By the end of the simulation, 
the torque experienced by the protoplanet will be positive 
on average, corresponding to outward migration on a time 
scale of ~ 2.8 Myr. 

If we assume that there exists a well defined long term 
mean for the torque experienced by the protoplanet, then 
we can express the torque as the sum of this mean torque 
and a fluctuating component: 



T{t)^T + Tj{t) 



(3) 



where T{t) is the torque experienced at time t, T is the mean 
torque, and Tf{t) is the fluctuating component. The running 
time average of the torque is given by 



IoT{t)dt 



then we can write 



{T)t 



T 



Tfit)dt. 



(4) 



(5) 



Assuming that the amplitude of the torque fluctuations 
Tf{t) about the mean T follow a Gaussian distribution 
with standard deviation ar and recur on a characteristic 
timescale tf, we can estimate from equation 5 that the fluc- 
tuations in the running mean satisfy 



\{Tt] 



(6) 



which may be used to estimate how long we need to run a 
simulation before we can expect the running time average of 
the torque to converge to a well defined value (i.e. when the 
fiuctuation estimate from equation 6 becomes significantly 
less than [r[). Inspection of figure 7 suggests that the tf 



Figure 9. This figure shows midplane density contours for the 
run G2. Note that the presence of the protoplanet is just de- 
tectable although the perturbations generated by the turbulence 
are of higher amplitude than those generated by the protoplanet. 
The protoplanet is located at {xp,yp)= (-3,0). 



is is typically about one half of a planetary orbit, Pp/2. 
Equation 6 tells us that the larger the relative amplitude of 
fiuctuations, the longer we have to integrate for the running 
mean to converge. 

It is clear that the level of fiuctuations associated with 
the torques is very much larger than the running mean in 
figure 7. If we take the mean value of the torque to be T ~ 
9 X 10~^ from the straight line in fig 7, then the typical 
magnitude of the fluctuations can be estimated to be ctt — 
2 X 10""", but with extreme fluctuations being a factor of 
~ 40 times larger than the mean during certain periods of 
the evolution. If we set or — 2 x 10~^, then the expected 
time for convergence is approximately 250 planetary orbits. 
This time scale is much longer than the time over which 
it is feasible to run simulations of this kind at the present 
time, and may provide at least a partial explanation of why 
convergence is not obtained in figure 8. But note that it 
may give a significant underestimate for the time required 
for convergence if the statistics of the fiuctuations are more 
complex involving say a superposition of fluctuations with 
varying characteristic times up to the total evolution time 
of the global disc. In such a situation convergence to small 
means may be practically impossible. 
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Figure 10. This figure shows midplane density contours for the run G2 in the region close to the protoplanet. The protoplanet is located 
at {rp,4>p)= (3,7r). The panels correspond to times 429.0687, 496.8698, 542.0710, and 634.4136, respectively. Note that the presence of 
the protoplanet is just detectable, although the perturbation it makes to the disc is of lower amplitude than the perturbations generated 
by the turbulence. 
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3.3 Global Model G2 

The protoplanet mass in this case was equivalent to ~ 10 
Earth masses. A snap shot of the midplane density distribu- 
tion generated by run G2 is shown in figure 9. The presence 
of the protoplanet is just discernible at {x, y)= (-3.,0). Close 
up images of the density field in the vicinity of the proto- 
planet are shown in figure 10, and again the presence of the 
protoplanet can just be detected along with the spiral waves 
excited by the protoplanet. The time variability of the den- 
sity field close to the protoplanet is apparent in this figure. 
Figure 11 shows the density field in the vicinity of the proto- 
planet in a laminar disc run that is otherwise equivalent to 
run G2, and it is worth noting that the density field shown 
in figure 11 remains essentially time independent once es- 
tablished. 

Figure 12 shows the torque per unit mass as a function 
of time for a laminar run equivalent to G2. The approximate 
constancy of the middle line in this plot (which is the total 
torque per unit mass on the protoplanet) shows that a well 
defined torque and migration time can be ascribed in this 
case. Figure 13 shows an equivalent plot for run G2. As in 
run Gl, the torque evolution in this case is characterised by 
large fluctuations about a small mean, such that the orbital 
evolution of the protoplanet is likely to occur as a random 
walk. 

The running time average of the torques for run G2 are 
plotted in figure 14. The time average of the total torque 
is again found to not converge, showing variations that in- 
dicate inward migration on average after ~ 7 orbits have 
elapsed, outward migration on average after ~ 18 orbits, 
and inward migration on average again at the very end of 
the simulation. Using similar arguments presented in sec- 
tion 3.2 to estimate the time for the running mean to con- 
verge, we again estimate that the run time required will be 
around 70 - 80 planet orbits. This is significantly more than 
the run time over which we are able to compute these mod- 
els at the present time. However, the torques from the two 
sides of the disc exhibit somewhat relatively smaller ampli- 
tude fluctuations when compared to simulation Gl. This is 
reflected in the shorter time estimated for convergence. It is 
also part of a tendency to have smaller relative fluctuations 
for larger protoplanet masses which makes net torque mea- 
surement easier. The same trend is seen in the shearing box 
simulations (see below). We also comment that the typical 
magnitude of the average torques from each side of the disc 
corresponds to the laminar value even in the low mass cases 
in spite of the large noise levels. 

Although that may be the case, large noise levels which 
apparently can occur on a variety of time scales will cause 
the migration of a low mass protoplanet embedded in a 
turbulent disc to differ substantially from that found in a 
laminar disc model, with periods of inward migration inter- 
spersed with outward migration. At this present time we are 
unable to give a reliable estimate of the net migration time 
of protoplanets in such disc models, or because there may be 
significant variations occurring on the longest evolutionary 
time scales of the global disc, give even a reliable indication 
of the direction of migration in the long term. Feedback of 
the orbital evolution on the turbulence may be significant in 
some cases. Numerical experiments are currently being per- 
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Figure 11. This figure shows midplane density contours for a 
laminar disc run with planet parameters identical to run G2. This 
figure should be compared with figure 10. 



formed to examine this and will be the subject of a future 
publication. 

In a recent paper, Winters, Balbus, & Hawley (2003a) 
have emphasised the 'chaotic' nature of MHD turbulence, 
with the implication being that the evolution of the tur- 
bulent state may be modified by small perturbations. The 
presence of a protoplanet induces such perturbations. Both 
simulations Gl and G2 were initiated with identical disc 
models (but not run G3), as commented upon in paper III. 
In paper III we noted that the global magnetic energy of 
the system was found to diverge (by a modest amount) once 
the planets had been inserted in the disc models. By com- 
paring figure 7 with 13, and 8 with 14, we can see the affect 
that the differing planet masses have on the local statistics 
of the turbulence. While the overall magnitude and form of 
the torque fluctuations are similar, the details are clearly 
different. Given that these fluctuations appear to strongly 
influence the orbital evolution of the protoplanets, we can 
conclude that the detailed orbital evolution of a low mass 
protoplanet in a turbulent disc is essentially indeterminate 
due to the feedback between the turbulence and the plane- 
tary perturbations. Interestingly enough though, we can say 
that even if the mean average torque results in inward mi- 
gration within a corresponding expected lifetime in the disc, 
large random fluctuations will result in some increased sur- 
vival probability for such a time for a subset of an ensemble 
of protoplanets. 
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Figure 12. This figure sliows tlie torque per unit mass exerted by 
the disc on the protoplanet for a laminar disc and planet param- 
eters identical to run G2. The upper line shows the torque due to 
the inner disc, the lower line shows the torque due to the outer 
disc, and the middle line shows the total torque. It is clear that 
a well defined torque is produced, with an associated migration 
time scale. 



Time Averaged Torques Per Unit Mass (Mp=10) 
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Figure 14. This figure shows the running time average of the 
torque per unit mass exerted by the disc for run G2. The upper 
line is the running time average of the torque acting on the planet 
due to the inner disc. The lower line is that due to the outer disc. 
The middle (not straight) line is the running time average of the 
total torque. The straight line is the total torque exerted on the 
protoplanet in a laminar disc run. We note that the total time 
averaged torque does not converge to a well defined value. 
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Figure 13. This figure shows the torque per unit mass exerted by 
the disc on the protoplanet for the run G2. The turbulence in this 
case generates strong fluctuations in the torque from each side of 
the disc, such that it becomes difficult to distinguish the torques 
arising from each side. The total torque fluctuates between pos- 
itive and negative values such that the associated migration will 
undergo a 'random walk'. 



3.4 Global Model G3 

Figure 15 shows a snap shot of the midplane density struc- 
ture in model G3, for which the protoplanet mass is equiv- 
alent to 30 Earth masses. The protoplanet is located at 
(a:,j/) = (-3.,0), and is clearly visible, along with the spiral 
waves it has generated. Figure 16 shows a series of close-up 
images of the midplane density near the protoplanet. Al- 
though the protoplanet and spiral waves it generates are 
clearly visible, it is also apparent that the density structure 
near the planet is highly variable. A close-up image of the 
density near the protoplanet for an equivalent laminar disc 



Mp=30 Earth Masses 



o 



Figure 15. This figure shows midplane density contours for the 
run G3. The presence of the protoplanet is clearly detectable, with 
the perturbations generated by the protoplanet being of similar 
magnitude to those generated by the turbulence. The protoplanet 
is located at {xp,yp)= (-3,0). 
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Figure 16. This figure shows midplane density contours for the run G3 in the region close to the protoplanet. The protoplanet is located 

at {rp,4>p)= (3,7r). The panels correspond to times 871.4805, 903.0731, 948.2260, and 1077.837, respectively. Note that the presence of the 

protoplanet is clearly detectable, with the perturbation it makes to the disc being of similar amplitude to the perturbations generated 

by the turbulence. Note that the appearance of the spiral waves generated by the protoplanet are time variable since they arc perturbed 
/ ^, , , , ^ ^ ^ ^ ^ (g) 0000 KAS, MNRAS OOCt; 000-000 

by the turbulence. ^ ' ' 



MHD turbulence and protoplanets 13 



Laminar Torques Per Unit Mass Vs Time (Mp = 30) 



Laminar Disc - Mp=30 Earth Masses 




1 xio 



-IxlO 



-2x10 



-,-5 
5 



-5 

-5 
-,-5 



-4x10 

5 

Time (Orbits) 

Figure 18. This figure sliows tlie torque per unit mass exerted by 
the disc on the protoplanet for a laminar disc and planet param- 
eters identical to run G3. The upper line shows the torque due to 
the inner disc, the lower line shows the torque due to the outer 
disc, and the middle line shows the total torque. It is clear that 
a well defined torque is produced, with an associated migration 
time scale. 
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Figure 17. This figure shows midplane density contours for a 
laminar disc run with planet parameters identical to run G3. This 
figure should be compared with figure 16. 



run is shown in figure 17 for comparison. Tlie density field 
and spiral wakes in tiiis case are essentially time indepen- 
dent. 

Tlie torque per unit mass as a function of time is shown 
in figure 18 for the laminar disc run equivalent to run G3. 
The middle line shows the net torque per unit mass, whose 
approximate constancy leads to a well defined migration 
time. The time evolution of the torque per unit mass for 
run G3 is shown in figure 19. A similar situation to that 
described for runs Gl and G2 is observed, with the torques 
being significantly modified by strong fiuctuations. Again 
the total torque can be seen to show periods of being pos- 
itive and negative, suggesting that the protoplanet would 
undergo a random walk rather than monotonic, inward mi- 
gration. 

We note, however, that as the protoplanet mass in- 
creases, and the spiral waves excited increase in amplitude 
to become comparable or larger that the turbulent density 
wakes, the torques due to the inner and outer disc begin 
to separate. This can be observed by comparing figures 7, 
13, and 19, which show the torques from the inner and outer 
discs becoming progressively distinguishable from each other 
as the planet mass increases. This is also accompanied by a 
reduction in the relative torque fluctuation noise level. 

Figure 20 shows the running time average of the torques 
per unit mass for run G3. The straight line in this figure 
is the time averaged total torque per unit mass obtained 
in the equivalent laminar disc run. As in runs Gl and G2, 
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Figure 19. This figure shows the torque per unit mass exerted by 
the disc on the protoplanet for the run G3. The turbulence in this 
case generates strong fiuctuations in the torque from each side of 
the disc, such that it becomes difficult to distinguish the torques 
arising from each side. The total torque fluctuates between pos- 
itive and negative values such that the associated migration will 
undergo a 'random walk'. 



the running time average of the total torque per unit mass 
fails to converge for the run time considered here. However, 
the relative fluctuation amplitudes are smaller in this case 
leading to a smaller anticipated time for convergence. This 
is part of the trend for larger mass protoplanets to produce 
larger amplitude perturbations that are more difficult for 
the turbulence to affect. 

If we consider the running mean of the total torque in 
figure 20 we can reasonably take a value of T ~ —6 x 10~®. 
A by-eye inspection of figure 19 indicates that the level of 
ffuctuation in the torques is ctt — 4 x 10~^. The predicted 
run time for convergence of the running mean from equa- 
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Figure 20. This figure shows the running time average of the 
torque per unit mass exerted by the disc for run G3. The upper 
line is the running time average of the torque acting on the planet 
due to the inner disc. The lower line is that due to the outer disc. 
The middle (not straight) line is the running time average of the 
total torque. The straight line is the total torque exerted on the 
protoplanet in a laminar disc run. We note that the total time 
averaged torque does not converge to a well defined value. 

tion 6 is then ~ 20 - 30 planetary orbits, which is shorter 
than for Gl or G2 and is comparable to the time for which 
the simulation has been run. 

Although it always corresponds to inward migration, 
the non convergence of the running mean toward a well de- 
fined value suggests that a simple picture of local turbulence 
in the vicinity of the planet, in which the state variables 
have well defined mean values on top of which are super- 
posed fluctuations with a well defined Gaussian spectrum, 
is not accurate. Instead, it appears that the global nature of 
the disc plays an important role in continuously modifying 
the local structure of the disc and turbulence. Communica- 
tion between different regions of the disc can affect the local 
properties such as the mean density, the amplitude, spatial, 
and temporal distribution of density and velocity fluctua- 
tions over long time scales, such that a local running mean 
is problematic to define. An examination of the global prop- 
erties of a turbulent disc model, such as the global magnetic 
energy (see e.g. papers I-III) show that there are short and 
longer time scale variations in the turbulence that refiect 
modifications to the local and global turbulence, and by 
implication the local disc structure. The results presented 
in section 3.6 for the local shearing box simulations show 
smaller relative fluctuations and greater convergence of the 
mean torques toward the expected value, indicating that the 
global properties of turbulent discs may play an important 
role by inducing longer time scale modification to the local 
state of the disc in the vicinity of an embedded protoplanet. 

3.5 Global Run G5 

The planet mass in global run G5 was equivalent to 3 Jupiter 
masses. The computational domain in this case was re- 
stricted so that the azimuthal interval ran between [0,7r/2] 
(thus allowing a reasonable run time), with the protoplanet 
being placed on a fixed circular orbit at (rp, 0p) = (2.5, 7i"/4). 



The physical parameters were identical to run G4, for which 
the azimuthal domain covered the full 2n. However, run G4 
was only run for ~ 11 planetary orbits, which is too short a 
time for anything other than the tendency for gap formation 
and the global disc morpholgy to be inferred. Consequently 
run G4 will not be discussed further here. 

Run G5 has been described in paper III, and showed 
a tendency toward clear gap formation with the response 
of the disc due to the presence of the planet being strongly 
non linear. Consequently the perturbations induced in the 
disc by the protoplanet are very much larger than those that 
arise because of the turbulence. This results in the fiuctu- 
ations in the torque experienced by the protoplanet being 
significantly smaller (in relative terms) than observed in the 
previously described runs Gl, G2, and G3, and a well defined 
running time average of the torque being obtained. Figure 21 
shows the running time average of the torque per unit mass 
obtained from run G5, with the upper line corresponding to 
the inner disc torque, the lowest line corresponding to the 
outer disc torque, and the middle line the running time av- 
erage of the total torque. It is clear from this figure that a 
large torque is exerted on the protoplanet prior to gap for- 
mation, but that as the gap proceeds to open and material 
is pushed away from the planet the torque diminishes. The 
running time averaged torques due to the inner and outer 
disc appear to be approaching well defined asymptotic val- 
ues, which are unaffected by turbulent fluctuations, but the 
continued decrease in the running time averages indicates 
that gap formation is still ongoing at the end of the simula- 
tion. 

We note that the use of a closed inner boundary in this 
simulation, combined with the close proximity of the planet 
to the inner boundary, cause the density of the inner disc to 
be maintained at an artificially high level after gap forma- 
tion. This leads to the net torque on the planet being neg- 
ative but close to zero. Under more general circumstances 
in which the inner disc can accrete onto the central star, 
an inner cavity is expected to form such that the torque on 
the protoplanet is dominated by the outer disc (e.g. Nel- 
son et al. 2000). If we adopt the disc model described in 
section 3.1 used to normalised the results presented in fig- 
ure 2, and estimate the migration time using equation 2 and 
a torque per unit mass due to the outer disc in figure 21 of 
T = — 10~^, then we obtain Tmig — 4 x lO** yr, for a planet 
at 5.2 AU. The type II migration time appropriate to gap 
forming protoplanets is given by the viscous evolution time 
Tmig = (2r-p)/(3i') whcre u — aH^Q, is the kinematic viscos- 
ity. For a disc model with a ~ 7 x 10"^ and H/r = 0.07, 
the estimated type II migration time is Tmig = 4.5 x 10* yr, 
in reasonable agreement with the result obtained from the 
simulation G5. 

We note that the above estimates for type II migration 
times correspond to disc models with full 2-n azimuthal do- 
mains, and that it is unclear which precise value the running 
mean of the outer disc torque will approach once gap for- 
mation is complete. Nonetheless, the reasonable agreement 
obtained in the estimates suggests that gap forming proto- 
planets in turbulent discs undergo migration at the expected 
type II rate. A similar result was obtained in paper II for 
type II migration rates in turbulent discs with full 2-n az- 
imuthal domains. 

It is clear that a well defined trend arises when consid- 
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Figure 21. This figure shows the running time average of the 
torque per unit mass exerted by the disc for run G5. The upper 
line is the running time average of the torque acting on the planet 
due to the inner disc. The lower line is that due to the outer 
disc. The middle line is the running time average of the total 
torque. The running time average for both the outer and inner 
torque converge to well defined values in this run, for which the 
interaction is non linear and the planet opens up a large gap. We 
note that our use of a closed inner boundary causes the density 
of the inner disc to be artificially increased after gap formation, 
such that the torque due to the inner disc is probably an over 
estimate, leading to a larger torque contribution. 

ering the interaction between embedded protoplanets and 
turbulent discs. Lower mass objects that are unable to per- 
turb the turbulent background flow significantly are subject 
to strong torque fluctuations that are likely to dominate 
their orbital evolution. As the protoplanet mass increases 
so that the amplitude of the spiral wakes that it excites 
become larger than the turbulent density fluctuations, the 
relative magnitudes of the torque fluctuations decrease, and 
the migration is likely to become similar to type I migration 
(although with a signiflcant noise component). For larger 
protoplanet masses that allow gap formation, the effect of 
the turbulent fluctuations is small, with the migration being 
essentially the same as the standard type II picture. These 
trends are also observed in the shearing box simulations that 
are described below. 

3.6 Shearing Box Simulation Results 

Details of the shearing box simulations Bal - Ba4 are 
given in table 2. These were each continued from a simu- 
lation with fully developed turbulence BaO after inserting 
a protoplanet with values of the dimensionless parameter 
GMp/{H^Q.l) = MpR^ /{M,H^) measuring the mass of the 
protoplanet equal to 0.1,0.3,1 and 2 respectively. In this 
section, Q.p and R axe the angular velocity and radius of 
the centre of the box. Thus simulations Bal and Ba2 are 
directly comparable to the global simulations G2 and G3 in 
terms of the relative strength of the protoplanet's perturba- 
tion on the disc. As in those cases the protoplanets remain 
embedded with no indication of gap formation. However, 
simulations BaS and Ba4 show gap formation and a non lin- 
ear perturbation of the disc as do simulations G4 and G5. 



We here discuss the force exerted on the protoplanet by the 
disc. In particular we study the component in the azimuthal 
direction Fy . Adopting a unit of length equal to the radius 
of the center of the box, this is also the torque. 

Density contour plots taken in the box mid plane near 
the end of the simulations Bal and Ba2 which have embed- 
ded protoplanets are shown in figure 22. As in the global 
simulations, the wakes produced by these protoplanets can 
be discerned even though the the medium is turbulent and 
produces erratic perturbations of them. Thus the tidal per- 
turbation of the disc produces measurable effects. 

Figure 23 shows running time averages of the torques 
acting on the protoplanet in simulations Bal-Ba4. Apart 
from simulation Bal, these commence from when the pro- 
toplanets were introduced into the simulations. In the case 
of Bal, the running average was started somewhat later. 

The averages were calculated over time periods of typ- 
ically 50 orbits at the centre of the box. In each case the 
average torques from the regions exterior to and interior to 
the protoplanet drag and accelerate the protoplanet as ex- 
pected. Although, as in the global simulations, fiuctuations 
in the one sided torques can be very large amounting to an 
order of magnitude or more greater than the typical aver- 
age value, the averages tend to approach reasonably steady 
values. Note that in a shearing box, symmetry considera- 
tions require that the mean torques from the two sides of the 
disc should ultimately be equal and opposite. However, some 
noise remains even after fifty orbits. The noise is stronger in 
the embedded cases and remains at the five to ten percent 
level after fifty orbits. The resulting mean net torque which 
would be zero in a laminar simulation suffers corresponding 
fluctuations. However, the effects of noise seems much less 
severe in simulations Ba3 and Ba4 in which the disc is non 
linearly perturbed, in agreement with the global run G5. 
In these cases turbulent fluctuations seem to be less able 
to modulate the torque estimates and non zero net torques 
arising from features added to bias the box are more readily 
measurable. 

To illustrate the reduced effects of noise on the simu- 
lations with strongly perturbing protoplanets, we compare 
running time averages of the torques acting on the proto- 
planet in simulation Ba4 with those obtained from the cor- 
responding laminar disc simulation in figure 24. The aver- 
ages commence from when the protoplanet was introduced. 
Corresponding plots from the two simulations are very sim- 
ilar. Note that the torques tend to relatively weaken at later 
times in the laminar case. This is because this simulation 
with {MpR^ / {Mt.H^) = 2 produces a wider and deeper gap 
than the turbulent case (see paper HI). This in turn results 
in less matter for the protoplanet to interact with locally on 
average and hence weaker average torques at later times. 

Finally we comment on the magnitudes of the time av- 
eraged one sided torques. As for the global runs these are 
similar for turbulent and corresponding laminar disc simu- 
lations. For the shearing box simulations, the fiducial value 
expected for the one sided linear torques in the 2D laminar 
case with no potential softening is 

Fy. = ^A.^R'nl[0,^^, (7) 

where the surface density E = {p)H, with (p) being the vol- 
ume and time averaged density in the box (see Ward 1997). 
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Figure 22. Density contour plots taken in the box mid plane for simulations Bal and Ba2 which have the smaller mass and hence fully 
embedded protoplanets. The wakes produced by these protoplanets can be discerned even though the medium is turbulent. 
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Figure 24. This figure shows running time averages of the 
torques acting on the protoplanet in simulation Ba4 commenc- 
ing from when the protoplanet was introduced. The uppermost 
curve extending to larger times gives the net torque acting on 
the planet from the inner regions while the corresponding low- 
ermost curve gives the contribution from the outer regions. The 
central short dashed curve gives the net torque. For comparison 
results from a laminar simulation with the same parameters but 
no magnetic field are presented. Corresponding plots from the 
two simulations can be almost superposed but the laminar case 
terminates at an earlier time. 

For simulation Bal, Fyo = 10~^^ in our computational units 
making the typical averaged torques about fifty percent of 
the fiducial value. Such a reduction may occur as a result 
of softening the gravitational potential due to the proto- 
planet (see section 3.1 and paper III). This would indicate 
that the magnitude of the torques is essentially given by the 
laminar disc theory. However, the turbulence adds a signif- 



icant noise component when the protoplanet mass is small 
enough to place the response in the linear regime. We com- 
ment that the averaged torque magnitudes in Bal and Ba2 
are consistent with a scaling oc Mp appropriate to the linear 
response regime but that the torque is weaker than sug- 
gested by such a scaling in Ba3 and Ba4. This is consistent 
with the response being non linear in those cases. 



4 DISCUSSION 

In this paper we have performed both global and local sim- 
ulations of embedded protoplanets interacting with a turbu- 
lent disc and studied the behaviour of the torques exerted 
between the disc and protoplanets and the consequences for 
orbital migration. The global simulations were for a disc 
with H/r = 0.07 that exhibited MHD turbulence with zero 
net magnetic flux with mean a ~ 0.007. The protoplanet 
masses considered were Mp = 3, 10 and 30 Earth masses, 
and 3 Jupiter masses. The local shearing box simulations 
can be characterized by values of the dimensionless param- 
eter MpR^/{MtH^). The simulations adopted 0.1,0.3,1.0, 
and 2.0 respectively. The first two of these are directly com- 
parable to the global simulations with Mp — lOM®, and 
A'lp — 3OM0 respectively. The latter two have gap forming 
protoplanets, but whose masses correspond to less massive 
planets than the 3 Jupiter mass planet considered in the 
global run G5, and enable the behaviour of the torques in 
the non linear gap forming regime to be studied. For a first 
study, the protoplanets considered here were held in fixed 
circular orbit. 

It was always found that the instantaneous torque expe- 
rienced by a protoplanet was a highly variable quantity on 
account of the protoplanet interacting with the turbulent 
density wakes that shear past it. For low mass protoplan- 
ets that are not able to begin to form a gap, the torque is 
dominated by these fluctuations, such that at any particu- 
lar time, the usual distinction between inner (positive) and 
outer (negative) disc torques is blurred. The net torque expe- 
rienced by embedded protoplanets oscillates between nega- 
tive and positive values, such that the protoplanet migration 
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Figure 23. This figure shows running time averages of the torques aeting on tlie protoplanet in simulations Bal top left panel, Ba2 top 
right panel, Ba,3 bottom left panel and Ba4 bottom rigfit panel. In each panel the upper lower and middle curves correspond to running 
averages of the torques acting on the planet due to the regions of the box interior to the planet, the torque due to the regions exterior 
to the planet, and the net torque respectively. 



is likely to occur as a random walk. This is in contrast to 
the monotonic inward drift normally associated with type I 
migration. 

In order to average out the erratic behaviour of the in- 
stantaneous torque a running time average is considered. We 
considered contributions from the inner and outer disc. We 
have been able to do this over a 20 orbital period timescale. 
Although these running averages took on values character- 
istically expected for type I migration, large fluctuations re- 
mained such that the net torques failed to converge to well 



defined values over the run times that arc currently feasi- 
ble. This is in contrast to laminar disc simulations for which 
convergence is achieved on an orbital time scale. 

Fluctuations in the averaged one sided and net torques 
were most severe for the smaller embedded masses. Thus 
for the 3M0 global simulation, fiuctuations in the one sided 
torques averaged over ~ 20 orbits were comparable to the 
averaged torques themselves resulting in an indeterminate 
direction of migration. However, in the case of the SOM® 
global simulation, fluctuations in the one sided averaged 
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torques were relatively smaller at about 20 percent of the 
mean values resulting in a more clearly delineated but still 
erratic inward migration. 

The existence of large fluctuations in the torques ex- 
erted by the disc on embedded low mass protoplanets and 
lack of knowledge of how these behave on long time scales 
makes definitive statements about the direction and rate of 
migration difficult to make. However, even in spite of large 
instantaneous fluctuations of more than one order of mag- 
nitude larger than the mean values, characteristic values of 
one sided torques averaged over a 20 orbit time scale are 
consistent with laminar type I estimates. Thus, when this 
estimate gives inward migration to the central object, we 
might expect the fluctuations to result in increased survival 
prospects for a subset of a statistical ensemble of embedded 
protoplanets. 

The behaviour of the fluctuations on long timescales re- 
mains an important issue that is impractical to explore at 
the present time. If there are significant fluctuations gener- 
ated in the global simulations that occur on the longest evo- 
lutionary time scales, convergence of torque running means 
becomes for practical purposes impossible to achieve and the 
migrational behaviour of low mass protoplanets considered 
as an ensemble would be very different from predictions of 
type I migration theory. 

We found that noise levels were relatively smaller in 
the local shearing box simulations. This gives some support 
to the idea of the existence of global fluctuations with long 
characteristic times in the global simulations which are ex- 
pected to be absent in the box simulations. 

We comment that the zero net magnetic flux models 
that we consider here generally give rise to turbulent discs 
that are more quiescent than those which arise when the 
magnetic flelds has net flux, especially if a net poloidal fleld 
is present (e.g. liawlcy, Gammic, & Balbus 1996; Hawley 
2001; Steinacker & Papaloizou 2002). Strong density fluc- 
tuations and contrasts, including persistent gaps, may arise 
in such discs, and would have a profound impact on the 
migration of embedded protoplanets. In particular, if MHD 
turbulence gives rise to global disc structures that lead to 
varying migration rates as a function of position in the disc, 
then migration of embedded protoplanets into these low- 
migration regions is likely to occur, increasing the lifetime 
of these planets in the disc. 

The results of both the global and local simulations 
show the same trend as a function of planet mass. For 
very low mass protoplanets the turbulent density wakes have 
higher amplitude than the spiral wake generated by the pro- 
toplanet. The torques exerted on the protoplanct are then 
such that the instantaneous turbulent fluctuations are very 
much larger than the running means. However, as the planet 
mass increases, the perturbation it makes on the disc starts 
to dominate its local neighbourhood. The expected sepa- 
ration between inner and outer disc torques becomes ap- 
parent and fluctuations become smaller relative to the run- 
ning mean torques. Eventually gap formation occurs. At this 
point we find a weakening of the torques exerted by the disc 
on the protoplanet due to the evacuation of gap material. 
There is then a transition to type II migration at a rate 
determined by the angular momentum transport in the dis- 
tant parts of the disc unaffected by the protoplanet. This 



has already been seen in the simulation of a 5 Jupiter mass 
protoplanet described in paper II. 
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